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Abstract 

We describe a PDE model of bedrock abrasion by impact of moving particles and show that by 
assuming unidirectional impacts the modification of a geometrical PDE due to Bloore Q exhibits 
circular arcs as solitary profiles. We demonstrate the existence and stability of these stationary, trav- 
elling shapes by numerical experiments based on finite difference approximations. Our simulations 
show that .depending on initial profile shape and other parameters, these circular profiles may evolve 
via long transients which, in a geological setting, may appear as non-circular stationary profiles. 

1 Introduction 

1.1 Motivation and references 

Bedrock abrasion by impact of moving particles is a dominant process shaping obstacles in river channels 
OX Abrasion is mostly due to saltating particles, a broad description of their mechanistic properties is 
presented in [0. Our goal is to identify and predict the geometry of obstacle profiles created by this 
process. In a recent paper, a discrete random model describing bedrock profile abrasion was developed 
y—{ ■ and its predictions compared with experiment [3J. The main observation was that stationary profiles 
emerged both in the laboratory experiments and the matching discrete simulations. This phenomenon 
has also been observed in Nature [4j,[5J,[9|. The aim of the present paper is to present a simple partial 
differential equation capturing the essence of the physical process and investigate whether it admits stable 
travelling front solutions which are compatible with the numerical and experimental results obtained in 
and observed in MIM- 



1.2 Basic assumptions and notations 

If (x, y, z) are Cartesian coordinates we take y along the vertical direction. We assume, for simplicity, 
that the system is uniform in the horizontal direction perpendicular the direction of motion of the abrading 
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Figure 1: Basic notations 



particles. Thus the bedrock profile starts off and remains independent of the z coordinate and is therefore 
given by a single function y = y(x, t) of space and time. The tangent to the profile makes an angle ip with 
respect to the horizontal, so the inward normal to the profile makes an angle if> — f with the horizontal. 
We also introduce 

y = dy/dx = tanip (1) 

and assume, following [3], that 

• small abraders are incident from the left making an angle <p with the horizontal. 

• the angle of the colliding face of large abraders is uniformly distributed between and |, i.e. the 
non-uniform flight of abraders does not change this term. 

As time proceeds, a point on the profile with constant x coordinate moves with vertical velocity given by 

_ dy - • ™ 

"^vertical — "777 ~~ Vi \ Z ) 

at 

the partial derivative being taken at fixed x. The velocity of the same point along the inward normal is: 

v n = -ycosip. (3) 

Abrasion due to collisions with incoming abraders with isotropic, uniform distribution has been first 
derived by Bloore Q, for the detailed derivation of the planar case see also H. Under the assumption of 
unidirectional abraders, Bloore's PDE can be readily modified by the inclination factor to obtain 

v n = (v + bp) sinO - 4>), (4) 

where p = —y"(l + y /2 )~ 3 / 2 is the geometric curvature and v, b are positive constants. Using the latter 
expression as well as ©, we obtain 

— y cos ip = (v + bp) sin('0 — <fi) = (v + bp) (sin if) cos — cos ifj sin 0) (5) 
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Using (OQ), from © we obtain 

— y = (v + bp) (cos 4>y' — sin 0) . (6) 

Without restricting generality, we assume = (i.e. that the x axis is aligned with the flight direction of 
the abraders) to obtain 

-y = y'(v + bp) = y'(v + by" (I + y' 2 )^). (V) 

Note that L = b/v is the only independent length scale of the equation. Rather than descibe the profile 
in terms of y considered as a function of x and t we may regard a; as a function of y and t: x = x(y,t), 
so we obtain: 

x= (v + bp) = {v + bx"{l+x' 2 y^) (8) 

where now x = d t x(y, t) keeping y fixed and x' = d y (y, t), keeping t fixed. 

Our equation is a model of collision-based abrasion, and the first (constant) term corresponds to the 
abrasion caused by small abrading particles whereas the second (curvature) term is significant if the 
abrader is large (cf. [Q,[|3],[]6)). Based on this physical background, in forthcoming numerical simula- 
tions we will assume that on non-convex parts of the profile the curvature term vanishes (as large abraders 
are not likely to hit concave domains). Accordingly, the modified equation reads: 

— V = y'( v + bp), where p = p if p > and p = if p < 0. (9) 



2 Existence of travelling waves 

Note that © is first order in time and second order in space and translationally invariant in time. It is also 
invariant under translations in space, both vertical and horizontal. Being first order in time (|7) defines 
a "flow" on the infinite dimensional space of profiles such that if we are given a smooth initial profile 
y(x, 0) = yinitiai(^) sa y> then, at least for some finite time y(x, t) is uniquely determined. In what follows 
we shall investigate whether, for smooth initial Cauchy data such that 2/mitiai(^) is monotonic increasing 
and tends rapidly to constant values for large negative and positive values of x, that the solution y(x, t) 
ultimately tends to an exact travelling front solution. If it does, then we shall derive the form it must take. 
If we make the common travelling wave ansatz 

y = f(x- a) (10) 

and after substituting into © solve the resulting ordinary differential equation for f(s) we find that the 
solutions are all of the form of portions of circles or straight lines which move to the right with speed 

c = v + - (11) 

a 

where a is the radius of the circle. Two special cases should be mentioned: For straight lines of course 
we have a = oo and c = v . In case of b = we have 

9y = _ v dy_ 
dt dx 

whose general solution is 

y = f(x — ct) , with c = v (13) 
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Figure 2: Circular arc with radius a = 2 and horizontal tangent at y — 1 evolving under © with 
L = b/v = 0.6. Numerical simulation by standard finite difference scheme, using N s = 300 spatial 
subdivisions and At = 10~ 5 timestep, profiles plotted after N t = 10 4 timesteps each. Observe that the 
profile moves with uniform translation, circular arc remains invariant under ©. 

thus if b = then any initial shape moves with to the right with speed c — v. These (local) travelling 
wave solutions travelling with speed given in (fTTI) may easily be derived from the alternative form of the 
governing equation © by making the ansatz x(y, t) = ct + x(y). 

The local solutions, made up of circles and straight lines may be patched together to give global 
solutions satisfying the boundary conditions. Matching conditions should be observed carefully. E.g. 
only segments with identical curvature may be patched together, otherwise the propagation speed c will 
suffer a jump. If we use boundary conditions representing a 'cliff (cf Figured) then we have two options: 

• (a) We satisfy horizontal tangency as boundary condition at both lower and upper end. In this case 
if we apply © prescribing that the curvature term vanishes on concave parts. In this case there is 
no exact travelling wave solutions as the two ends of the cliff can not be joined by a curve with 
constant curvature. The only exact solutions are half circles, however these contours are not a 
feasible representation of real cliffs. 

• (b) We satisfy the horizontal tangency only at the upper end. At the lower end we do not prescribe 
a boundary condition, i.e. we essentially cut off the profile at y — 0. This admits a one-parameter 
(R) family circular arcs which would be exact travelling solutions. This approach could be justified 
by considering that a vertex on the upper end is convex (thus it will be abraded by large abraders) 
whereas a vertex on the lower end is concave so large abraders to not smear it out. Small abraders, 
corresponding to the constant (Eikonal ) term do not abrade any vertices's. So, admitting a con- 
cave vertex at the lower end does not seem to contradict the physical assumptions. This type of 
boundary condition has been implemented numerically. Figure |2] shows the evolution of a circular 
arc under © and boundary condition (b). The numerical simulation was based on a standard finite 
difference scheme and we can observe that the circular arc remains invariant and moves by uniform 
translation. 
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3 Stability 



Next we investigate the stability of travelling waves. After finding analytical evidence in subsection 13.11 
for their local stability we run numerical simulation in sub section 13 .21 to test their global attractivity. 



3.1 Linear stability: analytical results 

Suppose y(x, t) is a solution of © when we assume the profile is convex, i.e. is negative, so that 

§ = , )£ (14) 



with b > 0. We consider a nearby solution y(x, t) + e(x, t) and substitute into (fl4)) and expand to lowest 
order in e we get 



(15) 



This equation is of the form 
with convection velocity 

and diffusivity 



de . . de , . d 2 e 



+ u(x,t)— = «0M)-^2 



u (x, t) = u - 6^ ^f- (17) 

K (x,t) = b 1 ± s-. (18) 

Essentially, the first term convects the disturbances with a positive space and time dependent speed u(x, t) 
and the second term is a diffusion term which causes the disturbance to spread out and decrease in 
amplitude. This is independent of the particular profile y(x, t). In case of straight lines 

y = A(x- vt) = As, (19) 

based on (fl7l) and (TT8l we find the convection velocity and diffusivity 

u(x,t)=v, K,(x,t) = & - - — . (20) 

We can aslo rewrite (fTol using or using the co-moving variables s, t 

de A d 2 e 

dt ~ 1 + A 2 ds 2 ( } 

which is the standard diffusion equation. So we can conclude that all travelling profiles are locally stable, 
moreover, for straight lines we have explicitly determined the diffusion and convection terms. 
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Figure 3: Initial profile given in (|22)) with A = 5 evolving under $7} with L = b/v = 0.6. Numerical 
simulation by standard finite difference scheme, using N s = 135 spatial subdivisions and At = 1CT 5 
timestep, profiles plotted after N t = 2 * 10 4 timesteps each (Since the shape changes faster initially, 
between the first an second profile an intermediate stage at 10.000 steps was inserted). Observe that 
profile approaches circular arc, next to last profile exact circle is plotted for better comparison. 

3.2 Nonlinear stability: numerical results 

Circular profiles are of key interest from the point of view of geological applications, we investigate their 
stability numerically. Figure |2] illustrated the local stability of circular arcs: by taking such an initial 
condition the circualr shape was preserved by the PDE. In te next step we consider the convex profile 



which satisfy y(0) = 0, y(oo) = 1 and we evolve these profile numerically considering the boundary 
condition (b) described in Section|2] The results are shown in Figure|3]and we can observe how the initial 
profile approaches a circular arc. For better comparison the last profile is plotted next to an exact circle. 

4 Comparison to experimental data 

Here we rely on laboratory experimental data obtained by Wilson [4J [8] using an annular, recirculating 
flume [8J. In an experiment with this flume [4], single cuboid marble blocks, 10,0 cm tall and 20 cm long, 
and spanning the 26 cm width of the flume, were fixed on the channel base as a flow obstruction. In our 
model this would correspond to a Heaviside step function as initial profile shape. The flume was loaded 
with sorted limestone pebbles to act as bedload abrasion tools, while water discharge was held constant 
to produce a flow speed of 3ms- 1. Under these flow conditions, the limestone pebbles moved by rolling 
and low saltations in a layer of one to two grain diameters depth along the flume base. Particles moved 
up the stoss (upstream facing) surface of an obstacle in one or several short hops, and launched off this 
face clearing the remainder of the obstacle. Resultant erosion of bedrock obstacles was measured using 
three dimensional laser scanning at intervals throughout experiment runs that lasted 400-690 minutes. 
Abrasion of obstacles was dominantly on the stoss surface whose initially square cross section evolved 
to an unpstream facing, almost convex surface. The lower 18.5 mm of the stoss surface were prone to 
edge chipping, and have been eliminated from consideration, the plotted profile is 81.5mm high. Rates 
of horizontal erosion were peaking near the top of the obstruction, and decreasing systematically towards 
the base, resulting in progressive reclining and convex-up rounding of the stoss side. Over time horizontal 
erosion converged to a common rate at all heights above the flume base, thus, the obstacle stoss faces 



Z/initiaiO) = tanh(Ax), x > 



(22) 
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Figure 4: Comparison of experimental and numerical data. Experimental data is based entirely on (01). 
Initial profile given in (l23l) evolving under © with L = b/v = 2.0. Numerical simulation by standard 
finite difference scheme, using iV s = 750 spatial subdivisions and At = 10 -6 timestep, profiles plotted 
after N t = 12500 timesteps each. Observe that numerical and experimental profile show very similar 
evolution. Profile approaches circular steady-state circular geometry with radius a = 131.5 only much 
beyond the range of the physical experiment. Location of inflection point is plotted both on experimental 
and numerical profiles. 

achieved a time-independent form that advanced downstream into the obstacle. These attributes are 
shared with bedrock obstacles in the Liwu River, Taiwan, that have been studied extensively IfTOl . We 
have set up numerical simulation applying our model to match the flume experiment. We approximated 
the initial profile by 

y = 41.5tanh(0.4z)+40 (23) 

The evolution of (|23l) under © is plotted in Figure ©a. Since this is a non-convex profile, we used the 
modified abrasion law defined in ®. We can observe that in this case the circular profiles are reached 
only after long transients, i.e. the inflection point (trajectory plotted in Figure HJa) lingers for an extended 
period of time above the x axis. This may account for the geologial observation of non-convex 'steady- 
state profiles which correspond actually to the long transients before the circular shape settles in. In 
Figure Hfc experiments are shown and Figure HJc. 
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5 Conclusions 

In this paper we set up the simple PDE model © for abrasion under unidirectional (parallel) impacts, 
based on Bloore's model for isotropic abrasion. By using the travelling wave ansatz (flOl) we found 
circular travelling solutions. After establishing their local (linear) stability analytically, we tested their 
global attractivity by finite difference simulations and found that they appear to be globally stable. We 
compared our numerical results to laboratory experiments of Wilson [4 J and found good agreement. Our 
simulation showed that nonconvex profiles may approach their final, circular steady-state geometry via 
long transients as inflection points may linger above the x axis for extended periods of time. This fact 
may account for the geological observation of (almost) stationary, non-convex profiles in riverbeds. 
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